Rejection sampling for complex models facing limited data
Ernesto Carrella
1 Introduction
Calibration finds the most likely model parameters given the data. When data is limited and the model complex, while likelihood can still be mechanically computed, it offers little value. The likelihood of the best parameters will be only insignificantly higher than that of many other candidates. Under these circumstances, calibration will output parameters that are neither credible nor account fully for the remaining uncertainty.
We propose giving up on finding the best parameter vector and instead focus on finding all the “good enough” ones. In practice we can achieve this by using rejection sampling(see chapter 14 in Russell and Norvig 2009; see also “rejection filtering” in Hartig et al. 2011). At its core rejection sampling simply means running the model many times with random parameters and rejecting all runs that do not match the data we observed.
The standard approach is to transform the set of all accepted model runs into confidence intervals (posteriors) around each parameter. Our main focus however is not the parameter space the accepted runs indicate, but the accepted runs themselves. The set of accepted model runs represent combinations of parameters, state variables and histories that match the data. We can use these runs for policy analysis (including management strategy evaluation): project each into the future and identify common risks and outputs.
Rejection filtering is popular in ecology as an operationalization of pattern oriented modeling(Grimm et al. 2005). Rejection filtering is also at the basis of the “ensemble ecosystem” approach(Baker, Gordon, and Bode 2017) where we generate a large set of qualitative Lotka-Volterra models and reject all those that do not match a set of known responses to shocks.
In fisheries population assessment science, many fish stocks around the world are data limited, hard to assess and manage. Catch-MSY (Martell and Froese 2013) is a fisheries application of rejection filtering. While the authors introduce it as an MCMC approach, it is better described as rejection approximate Bayesian computation (ABC) as in Beaumont, Zhang, and Balding (2002; see also Sisson 2018). Catch-MSY randomizes two parameters, logistic growth \(r\) and carrying capacity \(K\), and uses a complete time series of catches as input. Its only filter is the guessed interval for today’s depletion \(\frac{B_t}{K}\). Catch-MSY discovers the posterior for \(r\) and \(K\) by sampling them at random and rejecting all combinations for which simulated biomass is not within the depletion interval.
We differ from Catch-MSY in terms of objective. Catch-MSY rejection aims at building posteriors around its two parameters \(r\) and \(K\) and from them measuring maximum sustainable yield. This measure can be extracted out of the model for policy. We instead want to save the runs that generated those posteriors and simulate on each of them different policy shocks to study their overall robustness and rank their effectiveness.
In summary our objective is to sample from the space of acceptable fisheries without ever needing to identify it. We want to build a sampler of fisheries that are congruent with the data, the priors and the model we have. We want to do this without a likelihood function.
We explain the method in section 2 and then provide two examples. In section 3 we rank and optimize policies in a simple difference equation bio-economic model. In section 4 we show, inter alia, how a low SPR can be a predictor of imminent long-term economic improvement when contextualized in an agent-based model.
2 Methods
2.1 Rejection sampling
Given a set of evidence measured in the real world \(e_1,\dots,e_k\) and a simulation model \(m\) that takes as input parameters \(\Theta=\theta_1,\dots,\theta_n\) to generate the synthetic evidence \(\hat e_1(\Theta),\dots, \hat e_n (\Theta)\). Rejection sampling proceeds as follows:
- Draw a set of model parameters \(\hat\theta_1,\dots,\hat \theta_n\) from their prior distribution \(\pi(\Theta)\)
- Run the model and generate set of synthetic evidence \(\hat e_1(\hat\Theta),\dots, \hat e_n (\hat\Theta)\).
- If at least one synthetic evidence is not “close enough” to its real counterpart, reject (filter away) this simulation
We run through this loop until our computational budget is exhausted and collect all the simulations that were not rejected.
We use the ensemble of all accepted simulations to perform two kinds of tasks. First, we can study the posterior distribution of some parameter or other state variable. Studying the evolution of the posterior as we add evidence (rejecting more runs) allows us to understand how each evidence shapes the posterior. Second, we can apply the same policies to each fishery in the ensemble to compare their effect and predict their overall success.
2.2 Value of evidence
Models are often complicated and opaque. Model calibration even more so. It is often hard to link which piece of information changed what parameter by how much. With rejection sampling however one can simply compare the set of accepted fisheries before and after a piece of evidence is used as a filter. Because rejection sampling does no searching, the value of evidence problem is a lookup table operation.
As an example observe figure 2.1. In this fictitious example we are interested in \(x\) ( a parameter or state variable in the model) and how it changes as evidence is used. In this example there are two pieces of evidence and the prior on \(x\) is normal. The posterior, that is the set of \(x\) for the runs that pass both evidence filters, is bimodal. However had we only used evidence one, the posterior would have been still normal. Using only evidence two however gives the posterior its bimodality. Here then it is evidence two that is important for the shape of the posterior of \(x\).
Figure 2.1: A simple example to understand how to look at the effect of evidence in rejection sampling. A shows the distribution of \(x\) for all simulations (black) and for the simulations that passed both filters (yellow). B shows the distribution of \(x\) for runs that are accepted when we look only at evidence 1. C shows the same using only evidence 2. It is clear that the shape of the posterior is due to evidence 2. Evidence 1 rejects about half of all the initial simulations in the prior but does not alter the overall posterior distribution of \(x\) the way evidence 2 does
2.3 Policy analysis
In rejection sampling we use a theoretical model to generate synthetic evidence to filter out unrealistic fisheries. We can use the same model on all the accepted fisheries to predict their future. If the model allows it we can also simulate policy shocks and more complicated management procedures to study their average effect.
Figure 2.2 shows a trivial example. We accept only models that land within a range at time \(t=50\). The model runs that are accepted are then run for another 50 time steps to predict the future given our constraint.
Figure 2.2: An example of prediction with rejection sampling. We randomly produce 500 AR with random parameters and trend and reject all the runs that by year 50 are not within 10 and 15. Graphically the time series must cross the yellow segment. The lines that are accepted, colored in blue, are then run for another 50 time steps to predict the behaviour of the accepted models into the future.
3 Application 1: Bio-economic Model
3.1 Description
In this section we build a simple conceptual model to show how to assess the status of the stock and estimate policy effects. We show that this is possible in spite of the model remaining unidentified in all of its key parameters.
We use a bio-economic model Smith (1969; Seijo, Defeo, and Salas 1998) where a single logistic stock is targeted by boats whose number changes yearly in proportion to total profits. It can be summarised by two difference equations:
\[
\begin{align}
\Delta_t f &= \phi [ pqB_t - c]f_t \tag{3.1} \\
\Delta_t B &= r B_t [ 1- \frac{B_t}{K}]-qfB_t \tag{3.2}
\end{align}
\]
In this model, \(B_t\) is the biomass at year \(t\), \(f_t\) is the effort (expressed as number of boats), \(p\) is the revenue per ton caught, \(K\) is the carrying capacity, \(q\) is catchability and \(\phi\) represent how quickly boats enter or exit the fishery as a proportion of the profits observed. Table 3.1 lists the priors. Notice in particular the high carrying capacity: we expect this to be potentially a large fishery. Because we are not sure about how old the fishery is, the number of years to simulate for each run is an additional random parameter \(T\).
We want to condition this model on three pieces of information. First, we know that landings in the current year are somewhere between 10,000t and 15,000t. Second, we know that the fishery is currently not making profits. Third, we know that CPUE was higher ten years ago than it is today. This information cannot be sufficient to estimate even a simple bio-economic model. It is enough however to make meaningful policy analysis.
| Parameter | Minimum | Maximum | Meaning |
|---|---|---|---|
| \(p\) | 30$ | 90$ | revenue per ton of fish |
| \(c\) | 15,000$ | 60,000$ | annual cost per boat |
| \(q\) | 0.0002 | 0.0008 | catchability |
| \(K\) | 1,750,000t | 7,000,000t | carrying capacity |
| \(\phi\) | 0.0000025 | 0.00005 | profit to boat gain |
| r | 0.1 | 0.6 | logistic growth |
| T | 40 | 65 | age of the fishery |
3.2 Value of evidence
We run the model 1,000,000 times1. Each run uses a random set of parameters drawn from table 3.1. We discard all runs that do not match our evidence. That is, we discard runs where catches at year \(T\) are not between 10,000 and 15,000t; we discard runs where profits at time \(T\) are positive and we discard runs where \(\text{CPUE}_T>\text{CPUE}_{T-10}\). We accept 3,256 runs.
Even though our acceptance rate is 0.3% this is still not enough to identify any of the parameters. Figure 3.1 shows how none of the parameters can be identified given the data we have. Even worse, figure 3.2 proves that no MSY can be identified from the joint distribution of the logistic parameters. This is how it should be given the data limitations: MSY cannot be identified from one year of landings and two qualitative indicators.
Figure 3.1: The marginal posterior distribution for each of the parameters. While the distributions are not necessarily uniform as the prior was, there is non-zero density through the entire domain of each of the parameters. There is therefore no unique set of parameters identified by the data.
Figure 3.2: A scatter-plot where each point is an accepted run. Visually the density approximates the joint posterior distribution for logistic \(r\) and \(K\). Again, no real identification is possible as fisheries have been accepted through the entire joint domain. It is therefore not possible to estimate MSY and cognates given the evidence at hand.
Even though no parameter can be identified, we still generate informative posteriors on some state variables, in particular this year depletion \(\frac{B_T}{K}\). In figure 3.3 we show the posterior distribution of depletion as we incrementally filter away simulations that do not match evidence. Total landings in tons alone do not help identifying the status of the stock. This is somewhat unexpected because landings are low compared to our high carrying capacity prior. However it is possible for landings to be low if we are barely scratching the surface of a large untapped stock. The domain of the depletion posterior still ranges from 0% to 100%. Once we add the second piece of information (profits this year are negative) the posterior of depletion shifts left: we know now that we are dealing with a depleted stock. Adding the third piece of information (CPUE today is lower than ten years ago) shifts the posterior further to the left and it is clear that we deal with a heavily depleted stock.
Figure 3.3: The distribution of this year depletion \(\frac{B_T}{K}\) for all accepted simulations. A shows the prior: all simulations are accepted without looking at evidence; No depletion pattern is evident. B shows all the simulations whose last year’s catch conforms to evidence; it is clear that this piece of information alone does not help with determining depletion. C shows the depletion of simulations after filtering for both catches and profits; we can now see that the fishery is overfished. D shows the depletion after filtering for catches, negative profits and CPUE being higher ten years ago; it is clear that the fishery is heavily depleted.
Negative profits and a decade-old fall in CPUE are critical pieces of information then. Because the model we are using is a simple logistic, the drop in CPUE informs us that biomass was higher ten years ago. Falling profits implies that the fishery is either struggling with over-capacity or undershooting due to an earlier boom (perhaps connected to the larger CPUE observed ten years ago). Combining these with high prior we had on \(K\) and the low posterior we observe in depletion \(\frac{B_T}{K}\) implies that the current catches are small compared to what is achievable if the fishery is managed correctly.
3.3 Policy analysis
We would like to know which policy works best given the information we have. Because the model is not identified we cannot use parameters to identify a priori which policy will work. Instead we run each policy on all 3,256 accepted fisheries.
We simulate each policy applied to each fishery for 50 simulated years and collect averages over all simulated years and all fisheries. In table 3.2 we report average landings, profits and the percentage of years with depletion above 25% and 50% (which corresponds \(B_{\text{MSY}}\)) for each simulated policy. Business as usual (no policy) lands large quantities of fish, but does so inefficiently as profits are overall negative. Most simulated fisheries remain trapped in a continuous boom-bust cycle. Setting a TAC as catches at year \(T\) (last year before policy) proves too restrictive: landings are minimal even after biomass recovers. This makes sense: we know the fishery is heavily depleted and current catches are likely far below average values. Setting a maximum effort to 80% of boats at time \(T\) provides only moderately more landings than business as usual but larger profits as the policy prevents over-capacity booms. Setting TAC to each fishery’s theoretical best (assuming therefore some form of miraculous identification) provides the most landings but not necessarily the most profits as the fishery can still become crowded with many boats chasing an ever smaller share of the theoretically sound TAC.
| Average policy results | ||||
|---|---|---|---|---|
| Total Landings | Total Profits | 0.25 Depletion | 0.5 Depletion | |
| Business as usual | 6,497,560t | −$5,292,318 | 30.43% | 13.59% |
| 80% current boats | 7,804,540t | $286,430,796 | 56.67% | 36.90% |
| Last year catch TAC | 432,365t | −$57,284,877 | 60.21% | 51.12% |
| MSY TAC | 9,561,784t | $68,281,386 | 57.93% | 43.21% |
| Closed | 0t | $0 | 69.78% | 59.99% |
Besides policy analysis we can manipulate the accepted fisheries as with any other simulation model to look for a policy that maximizes our objective(Bailey et al. 2018). Focusing on effort control we look for the number of boats that would maximize profits over the next 50 years. Knowing the exact number of boats that would maximize profits requires a fully identified model. We instead set maximum effort to be a percentage of boats we observe this year (time \(T\)) and look for maximum profits “on average” across all accepted fisheries.
In table 3.3 we show the result of a grid search looking for the right % of current boats to set as maximum effort. Profits are maximized when setting maximum effort to between 50% and 60% of current effort.
| Policy optimization | ||||
|---|---|---|---|---|
| Looking for the maximum effort to impose as % of boats currently in the system | ||||
| Total Landings | Total Profits | 0.5BMSY | BMSY | |
| 10.0% current boats | 2,359,611t | $122,542,205 | 68.65% | 58.42% |
| 20.0% current boats | 4,155,692t | $208,558,019 | 67.43% | 56.75% |
| 30.0% current boats | 5,468,769t | $263,642,642 | 66.17% | 54.20% |
| 40.0% current boats | 6,390,060t | $294,053,626 | 64.72% | 50.59% |
| 50.0% current boats | 7,009,285t | $306,049,585 | 62.89% | 46.78% |
| 60.0% current boats | 7,409,192t | $305,568,079 | 60.84% | 42.71% |
| 70.0% current boats | 7,656,195t | $297,696,871 | 58.70% | 39.45% |
| 80.0% current boats | 7,804,540t | $286,430,796 | 56.67% | 36.90% |
| 90.0% current boats | 7,914,020t | $275,320,579 | 54.74% | 34.99% |
| 100.0% current boats | 8,014,010t | $266,004,973 | 53.09% | 33.31% |
For complicated multi-parameter policies one could reduce overfitting and get more accurate estimates of policy effects by generating a second set of accepted fisheries to use as a testing set.
Profits and landings were used and optimized in absolute terms. This is because our first piece of evidence (catches between 10,000t and 15,000t) keeps the set of acceptable fisheries approximately comparable. If no such evidence were available landings and profits would be better normalized or studied as relative to a baseline in order to keep averages from being affected by outliers.
4 Application 2: hairtail fishing in the Bungo channel
4.1 Description
In this example we model the Japanese coastal trolling fishery based in Usuki as described in Makino et al. (2017). As of 2011 it was composed of 45 boats, all smaller than 5 gross tons, targeting hairtail (Trichiurus japonicus). Sale price of hairtail depends on its size and letting it grow for 6 months can triple its value. The contribution of the Usuki trolling fishery to the overall fishing mortality is however limited as trolling boats share the stock with purse seiners who target other species and land hairtail as a byproduct. The fishery is not data-limited as Watari et al. (2017) presents a full stock assessment. The purpose of this example is one of data-degradation: what changes when we pretend to have less data than what is actually available.
Makino et al. (2017) focuses on the economic side of the fishery: number of fishers is in decline, most of them earn about 5-6M JPY which covers fixed and living costs but is not enough to save the 25M JPY needed to build a new boat. The focus on individual boats makes it ideal for an agent-based model analysis and here we create a simple POSEIDON (Bailey et al. 2018) application for it.
We sketch here the key ideas of the model but leave the full details in the appendix.
The biological operating model is length-based, with 5cm length-bins using the fixed boxcar method (Goudriaan 1986). Recruitment is noisy and happens twice a year (spring and autumn) assuming a Beverton-Holt stock-recruitment relationship.
Geographically fish is spread out over a 5-by-5 grid, each cell a square 12km wide.
We model each trolling boat separately. Boats can fish any day of the year and choose the fishing location by trial and error through the explore-exploit-imitate algorithm (Carrella, Bailey, and Madsen 2019).
We track the profits made by each boat. Those that manage to accumulate more cash than their target savings (after accounting for living and fixed costs) spawn a successor (that is, create a new boat that joins the fishery). Those that fail to cover their living costs will eventually leave.
The purse seiner fleet is not agentized and is instead a simple constant fishing mortality \(F\) applied to the stock at the end of each year.
We run the model for at most 45 years, starting the trolling fishery ab novo at year 1. We condition the model on only six pieces of information. First, total landings (trolling and purse seiners combined) in the fishery have never been above 15,000t. This is evident from figure 2 in Watari et al. (2017). Second, current landings of the Usuki trolling fleet is between 250t and 1,850t. The lower bound is used as an historical target in Makino et al. (2017), the upper bound is the non-net fishing mortality computed for 2011 in Watari et al. (2017). The upper bound is probably too high as it includes other sources of fishing mortality. Third, current SPR is unsustainable but not catastrophically so and we will accept fisheries whose current SPR is between .10 and .25. Watari et al. (2017) cites an SPR for the fishery of .21 in 2011. Fourth, the Usuki trolling fishery today has less than 60 active boats. Fifth, the Usuki trolling fishery represents at most 30% of the total landings for the stock. Sixth, the fishery is at least 30 years old.
Because our model is length rather than age based, we compute SPR via a length-based approximation, assuming median life-history parameters for the species (using Fishlife from Thorson et al. 2017; rFishBase from Boettiger, Lang, and Wainwright 2012). In other words, when computing SPR we do not assume we know the correct life-history parameters within the simulation and as such we expect this filter to be quite noisy.
This POSEIDON application is simple and lacks many details, inter alia it assumes knife-edge maturity, a simple geography, no seasonal variation in ex-vessel prices and simple stochasticity for the bi-annual spawning pulses. In spite of its simplicity it still contains too many parameters (see table 3.1) to hope for identification given the little data we condition with.
| Variable | Distribution | Meaning | Source |
|---|---|---|---|
| catchability | \(U[10^{-5},10^{-9}]\) | % of biomass caught per hour spent fishing | - |
| \(S_1\) | \(U6.2,9.3]\) | Logistic selectivity, first parameter | (size limit at 25cm) |
| \(S_2\) | \(U[0.17,0.26]\) | Logistic selectivity, second parameter | (size limit at 25cm) |
| Cell with most biomass | \(U[1,5] \times U [1,5]\) | Cell in the 5x5 grid with the most biomass | - |
| Biomass smoothing | \(U[0.1,1]\) | Parameter defining the geographical dispersion of biomass | - |
| Virgin recruits \(R_0\) | \(U[37M,62M]\) | Number of annual recruits at virgin biomass | \(R_0\) that would generate 10,000t to 200,000t of virgin biomass |
| \(\Phi\) | \(U[0.27,0.31]\) | Ratio \(R_0\) to spawning stock biomass | Implied by knife-edge maturity from Fishlife median \(L_{\text{mat}}\) |
| Steepness | \(U[0.2,0.99]\) | Proportion of current recruits to virgin when 20% of biomass left | - |
| Initial depletion | \(U[0.7,1]\) | Proportion of biomass to virgin available as simulation starts | - |
| Life history parameters \(L_{\infty},K,L_{\text{mat}},M\) | Drawn from Fishlife | Parameters governing the speed of growth for all fish | Drawn from Fishlife multi-normal distribution |
| Allometric parameters \(a,b\) | Bootstrap sampled from rFishbase | Parameters converting fish length to weight | Sampled with replacement from Fishbase |
| Hourly variable cost (JPY/hr) | \(U[100,140]\) | Monetary cost per hour spent at sea (fuel,engine) | Makino et al. (2017) |
| Hourly effort cost (JPY/hr) | \(U[700,800]\) | Monetary cost per hour spent fishing (film, bait and boxing) | Makino et al. (2017) |
| Hold size (kg) | \(U[1000,3000]\) | Maximum amount of fish transportable by a boat | 5GT boats |
| Savings target (JPY) | \(U[20M,30M]\) | Money to accumulate to spawn additional boat | Makino et al. (2017) |
| Yearly Expenses | \(U[4.5M,7M]\) | Yearly expenditure (living costs, plus fixed costs) | Makino et al. (2017) and census |
| Daily probability of fishing | \(U[.75,.9]\) | Probability of going fishing each day of the year | Comparing boat numbers with total effort recorded |
| Yearly exogenous fishing mortality | \(U[.3,.8]\) | Mortality rate due to other fishers targeting the stock | - |
4.2 Value of evidence
We run the model 354,775 times and accept only 743 fisheries. The extremely low acceptance rate (0.2%) is explained by the unnecessarily wide priors, particularly in catchability, steepness and exogenous fishing mortality. Many draws from the prior are fisheries that are not economically viable or biologically sustainable for more than a few years. In terms of marginal distributions most posteriors have not shifted from their priors but many of the parameters that were independent are now correlated.
The joint distribution of \(M\) (natural mortality) and \(K\) (growth coefficient) is particularly affected by rejection sampling. The prior of \(\frac{M}{K}\), which we inherit from Fishlife, is bell-shaped and centered around 1.3. The posterior is still bell-shaped but shifted to the left (median of 1.09).
We show in figure 4.1 how evidence would shift this posterior further. In subfigure B we draw the two hypothetical posteriors had we accepted SPR above .3 or below .1 instead of the .1-.25 interval we are actually using. Assuming high SPR moves the \(\frac{M}{K}\) posterior left, while low SPR moves posterior right and makes it bimodal. This in essence means that low SPR could be a symptom of low \(\frac{M}{K}\) (a fish that naturally grows slower and dies quicker) rather than any true signal of stock health. Sensitivity of length based SPR to mispecifications of \(\frac{M}{K}\) is well understood (Hordyk et al. 2014, 2019; Froese et al. 2019b, 2019a) and here we are simply picking up the conflict between Fishlife’s wide priors and using its modal values to compute SPR. In other words, on its own, length-based SPR could be either a sign of overfishing or a sign of mispecified life history parameters.
Figure 4.1: Ratio between natural mortality M and growth rate K; subfigure A shows prior and posterior distribution of the \(\frac{M}{K}\), the prior is drawn from FishLife and the posterior represents the fisheries that passed all filters, the posterior has shifted mostly to the left ;subfigure B compares the \(\frac{M}{K}\) posterior had we only SPR above 30 (cyan) or below 10 (green); subfigure C shows the posterior achieved by accepting only fisheries who have more than 60 boats active today (orange).
A far more counterintuitive effect of evidence is the limit on the number of boats (no more than 60 boats in Usuki today) as shown in figure 4.1, subfigure C. If we accept only runs with more than 60 boats, the \(\frac{M}{K}\) posterior is bi-modal and shifts further to the right: we expect a higher \(\frac{M}{K}\). This is counter-intuitive because higher \(\frac{M}{K}\) means that the fish grows slower and is less likely to get to maturity. We would expect this stock to be able to sustain only a smaller fleet.
What we uncovered is a collider bias (a Monty Hall effect as in Burns and Wieth 2004). It is true that higher \(\frac{M}{K}\) implies a more vulnerable fishery: the higher \(\frac{M}{K}\) in the prior the fewer years it takes for a random fishery to collapse. High \(\frac{M}{K}\) also implies less economic viability: larger fish is paid more and the higher \(\frac{M}{K}\) gets the less likely the fish is to get large enough to command a price premium. This however is a priori knowledge: a posteriori we are not looking at the whole space of fisheries but only the specific space of fisheries that survived until today (at least \(250t\) landed in Usuki this year) and still representing only a portion of total landings.
Fisheries with high \(\frac{M}{K}\) that survive are quantity-driven rather than value-driven. If fish matures slowly and is less likely to survive, boats need to catch it when it is shorter. If short fish sells for less, boats need to catch more of it per unit of effort. Knowing both that \(\frac{M}{K}\) is high and that the fishery survived implies that other model parameters must change to make it so. This is done mostly by exploiting the wide priors on weight-length parameters \(a,b\). In the prior there is no correlation between \(\frac{M}{K}\) and weight-length parameters, in the posterior the correlation is 0.88: in an accepted fishery with \(\frac{M}{K}=1.2\) a 50cm fish weighs 200g, but at \(\frac{M}{K}=2\) it has to weigh 600g. The reverse wouldn’t work: high \(\frac{M}{K}\) with low weight would be economically unviable, low \(\frac{M}{K}\) with high weight would exceed landings limits (either in absolute tonnage or relative to the purse seine fleet). Therefore, high \(\frac{M}{K}\) in the posterior implies a quantity-driven fishery with a high number of boats.
This collider effect can work to our advantage. Even with data, \(\frac{M}{K}\) and in particular natural mortality \(M\) are hard to estimate and drive many stock assessment results (Mangel et al. 2013). Here however we have connected the ratio to weight-length parameters which can be estimated cheaply. The heavier the fish, the higher our expected \(\frac{M}{K}\) becomes.
The same collider effect explains an even more puzzling prediction: runs with lower SPR will “do better” (attract more boats) in the future, as shown in figure 4.2. As explained above, low SPR is associated with high \(\frac{M}{K}\) and high \(\frac{M}{K}\) is associated with quantity-driven fisheries with many boats. High participation fisheries are usually filtered away because of our 60 boats limit but participation today may be temporarily low due to recent recruitment failures (critically 10 to 5 years ago). In context then, a low SPR may indicate a high participation fishery that is currently under-strength. This kind of fishery will eventually snap back to its long-term high participation equilibrium.
Figure 4.2: Projected and past realizations of participation and recruitment for accepted fisheries as well as fisheries that would have been accepted had SPR evidence been different. Each line is a separate simulation, bolded lines are median values. The subset of runs with low SPR and also M/K ratio above two have been grouped separately (in black). Low SPR today implies higher participation tomorrow if associated with higher M/K. This is better explained by looking at recruitment (subfigure B): these fisheries in equilibrium have higher participation but are over-correcting to low recruitment pulses 5 to 10 years ago and are likely to rebound.
We can generate another counter-intuitive result by focusing on evidence we did not use. Makino et al. (2017) mentions how the number of Usuki fishers have decreased by 40% in the last 10 years. We did not use this figure because it is the total loss across all fisheries, not hairtail specific. In this paragraph however, we add this as another piece of evidence. Again rejection sampling helps us contextualize this piece of information: a large drop in boats may mean an overfished stock, a series of weak recruitments, or a too aggressive undershooting in participation that will self-correct as profits return. In figure 4.3 we show future participation splitting the accepted fisheries between those that suffered at least a 40% boat reduction in the last ten years and the fisheries that didn’t. Within 15 years there is no real difference in participation and in fact fisheries that declined rapidly tend to see higher participation in the short run.
Figure 4.3: Projected fishery participation for the next 15 years. Each line is a separate accepted fishery, bolded lines are yearly medians. Accepted fisheries were split into two groups: those that saw at least a 40% decline in boats within the last ten years and those that didn’t. Large participation decline is not a factor predicting further drops in the long run and to the contrary may mean larger participation in the short run.
4.3 Policy analysis
For the accepted fisheries, business as usual will not result in any long-term change in the number of fishers or any other important indicator. This we show in figure 4.4. It remains however a depleted fishery with biomass at approximately 15% of carrying capacity, justifying additional policies.
Figure 4.4: Projected participation and biomass for the accepted fisheries when no policy is implemented. Each line is a simulation,the bolded line is the median yearly value. While each individual run will exhibit noise and short term trends, the overall dynamic is one of long-term equilibrium around current values.
We assume here that we can only enforce policies on the Usuki trolling fleet and cannot affect fishing mortality outside of it. Because we constrained the Usuki fishery to represent at most 30% of the landings, we cannot expect regulations to achieve much rebuilding.
We can still make a difference for the economic health of this fleet in two ways. First, we could try marginally improving biomass in order to increase CPUE and therefore profitability. Second, we could try changing selectivity to increase the average size of fish caught and therefore earning more money for each ton of fish landed.
We compare three policies. The first is just closing down the fishery to new entrants and re-entry. The second combines the first policy with imposing a limit of 200 days a year for fishing. The third is to combine the first policy with mandating a new gear that catches larger fish (specifically, increases by 10% \(S_1\) and decreases by 10% \(S_2\) of the logistic selectivity).
While seasonal closures do achieve higher profits per trip and higher selectivity increases revenue per unit caught, these benefits are negated by the purse seine non-Usuki fishery which increases its relative importance in landings. Seasonal closures generate higher profits per trip but these fail to cover fixed costs. If the overall objective is either increased profits or higher biomass, banning entry/re-entry alone performs just as well as more complex policies. We show these results in figure 4.5.
Figure 4.5: Key indicators projections for three policies (no entry or re-entry, 200 days fishing limit, higher selectivity); each line is a separate fishery, bolded lines are the median yearly value. Better selectivity increases revenue per kg caught (subfigure D) and fishing limit increases profits per trip (subfigure B). In general however all policies cause more landings to occur outside the Usuki trolling fishery (subfigure E) and selectivity improvements do not improve overall savings any better than just closing to new entrants (subfigure F)
4.4 Comparison to data-rich results
Biomass in the stock assessment produced in Watari et al. (2017) is estimated at 4,896t in 2011. Here the mean current biomass is 3,477t with standard deviation of 1,806t. We are about one standard deviation away using only limited evidence.
The main difference between our work and the data-rich simulations in Makino et al. (2017) is that our results are far more optimistic in the status quo. In the original data-rich work the unmanaged fishery constantly deteriorates, both biologically and economically. Our projections instead show a long term equilibrium around current values.
The discrepancy can be explained by the underlying differences in the economic model. In our paper boats quit when earnings dry up and bank balances turn negative. This helps the fishery recover over time and avoid collapses. In Makino et al. (2017) however, fishing mortality is a free parameter and as such fishing pressure can continue unabated.
One contributing factor for discrepancies is also our choice for life-history parameters. We used “off-the-shelf” FishLife distributions as priors. These are wide and not necessarily centered near the real values from the stock assessment. For example, in its stock assessment Watari et al. (2017) assumes maximum age for hairtail at 5 (the maximum ever observed in the area). Using Quinn and Deriso (1999) thumb rule (Hordyk et al. 2014) this would imply a natural mortality of at least 0.9. Fishlife average is 0.39 with a standard deviation of 0.11.
5 Discussion
5.1 Rejection sampling and the law of decreasing credibility
Fitting complex models to limited data generates identification issues(see Canova and Sala 2009). Some parameters may not affect the likelihood function (under-identification) or affect it only proportionally with others (partial identification). The curvature of the likelihood function around each parameter may be small enough to make maximization difficult in practice (weak identification).
The common solution to lack of identification is to add assumptions to achieve it with the data we have(see for example Manski 2007). These assumptions often involve data we couldn’t collect. For example, to run a virtual population analysis (VPA) we require years of observations to track fish cohorts through time. However if we assume that the population is in equilibrium and that future cohorts will resemble current ones we can run VPA with a single snapshot (see pseudo-cohort adjustment in Sparre and Venema 1998; Mildenberger, Taylor, and Wolff 2017). Assumptions however are subject to the “law of decreasing credibility”(Manski 2003): the more we want to reduce uncertainty the less realistic our assumptions become.
For certain questions however reducing the uncertainty on the parameter posterior is unnecessary. If we are interested in the average behaviour of the model given some observed data we may not require point estimation for the parameters. As long as we can sample from the set of fisheries that match the data, we can just average out their projections as the central limit theorem will bound their expectation. In fisheries this approach is championed by the DLM toolkit (Carruthers et al. 2014; Hordyk and Carruthers 2018).
It should be noted however that simulation does not avoid the “law of decreasing credibility”. The model that generates candidate fisheries contains many assumptions of varying degrees of credibility. The advantage however is that the model assumptions tend to be about the dynamics of the system and theoretical causal mechanisms (the way things work) rather than about data and evidence we do not possess (the way things are). The cost of simulation is that we can only focus on average behaviour of the model which in some applications may be not interesting.
5.2 Theory as a substitute for data
Varian (1989) summarised the economics approach to under-identification with the phrase “theory as a substitute for data”. When there are not enough observations for an extrapolation, theory can help connect the evidence we do have with the parameters we are interested in. This both in the sense of telling us what parameter matters and what observations can help identify it.
Simulation and rejection sampling is a way to automatically connect theory and data in a model-agnostic way in exchange for a steep computational price. It is a computational substitute to drawing causal directed acyclic graphs connecting evidence and parameters through conditioning (Pearl 2017, 2014).
A strand of simulation literature, particularly in economics and finance focuses on stylized facts (Cramer and Trimborn 2019; Franke and Westerhoff 2016). These are broad strokes descriptions of a system which are used to validate a model and its parametrization. Rejection sampling is different from the stylized facts approach, even when fed stylized evidence as rejection sampling better accounts for uncertainty. With limited evidence the problem is not how to build a model that replicates stylized facts, rather the problem is to account for all the many ways in which that is possible. Once we are able to sample from this space we can look for common patterns and consequences.
In the first application for instance, we could have identified the fishery as depleted by intuition alone without running a model and rejection sampling. The value added rejection sampling provided was a full posterior distribution and therefore a probabilistic description of the uncertainty of parameters and state variables. We can generate this posterior because we are using theory to establish causal mechanisms linking parameters and dynamics to the evidence and observations.
The cost of leaning on theoretical causal mechanisms is that filtering is only as good as the model we are using.
We face a tradeoff: we would like the model to be simple so that it runs fast and produces a large ensemble of accepted fisheries, but also complicated enough that our causal mechanisms are realistic (should we add recruitment noise? technological change in catchability?).
We also need to add enough complexity in the model to produce the evidence we need to filter against. For example, if we observe catch-at-length then our model must contain age or length structure.
As computing becomes cheaper, this tradeoff will tilt towards more complicated models.
5.3 Connection to Approximate Bayesian Computation
Hartig et al. (2011) lists rejection filtering as a separate method from other likelihood-free inferential methods. However it is better to consider rejection filtering as a special case of approximate Bayesian computation (ABC) where summary statistics are binary (pass or fail). Zhang et al. (2017) links ABC to pattern oriented modeling (which is often used as synonym for rejection filtering).
There are two advantages in tying rejection filtering to the rich ABC literature. First, the convergence of ABC methods is well understood and guaranteed under mild assumptions (Barber, Voss, and Webster 2015).
Second, the ABC methodology provides computationally more efficient alternatives to simple rejection sampling. One popular technique replaces random sampling with either MCMC or SMC techniques to reduce the number of runs needed (Beaumont et al. 2009). This proved unnecessary here but would become important with more evidence or slower models (Grazzini, Richiardi, and Tsionas 2017 shows how these techniques can actually be slower than simple rejection ABC when applied to simple problems).
5.4 Conclusion
We presented a method to use models given limited data for policy analysis accounting for identification failures and the resulting uncertainty. We presented two examples to show how the method is agnostic both to the kind of model used (difference equation or agent-based model) and the kind of evidence brought to bear. This method requires no likelihood function and can be applied in any circumstance, if we are able to pay the computational cost.
More importantly however, we proposed an approach to data-limited fishery analysis that does not involve making assumptions in order to stretch the data into an indicator to use directly for policy. Rather we proposed sampling from the set of fisheries that match the evidence and then look for conclusions that are robust enough to our lack of identification.
Bibliography
Bailey, Richard M., Ernesto Carrella, Robert Axtell, Matthew G. Burgess, Reniel B. Cabral, Michael Drexler, Chris Dorsett, Jens Koed Madsen, Andreas Merkl, and Steven Saul. 2018. “A computational approach to managing coupled human–environmental systems: the POSEIDON model of ocean fisheries.” Sustainability Science, June, 1–17. https://doi.org/10.1007/s11625-018-0579-9.
Baker, Christopher M., Ascelin Gordon, and Michael Bode. 2017. “Ensemble ecosystem modeling for predicting ecosystem response to predator reintroduction.” Conservation Biology 31 (2): 376–84. https://doi.org/10.1111/cobi.12798.
Barber, Stuart, Jochen Voss, and Mark Webster. 2015. “The rate of convergence for approximate Bayesian computation.” Electronic Journal of Statistics 9: 80–105. https://doi.org/10.1214/15-EJS988.
Beaumont, Mark A, Jean-Marie Cornuet, Jean-Michel Marin, and Christian P Robert. 2009. “Adaptive approximate Bayesian computation.” Biometrika 96 (4): 983–90. https://doi.org/10.1093/biomet/asp052.
Beaumont, Mark A, Wenyang Zhang, and David J Balding. 2002. “Approximate Bayesian computation in population genetics.” Genetics 162 (4): 2025–35. https://doi.org/Genetics December 1, 2002 vol. 162 no. 4 2025-2035.
Boettiger, C., D. T. Lang, and P. C. Wainwright. 2012. “Rfishbase: Exploring, manipulating and visualizing FishBase data from R.” Journal of Fish Biology 81 (6): 2030–9. https://doi.org/10.1111/j.1095-8649.2012.03464.x.
Burns, Bruce D., and Mareike Wieth. 2004. “The collider principle in causal reasoning: Why the Monty Hall dilemma is so hard.” Journal of Experimental Psychology: General 133 (3): 434–49. https://doi.org/10.1037/0096-3445.133.3.434.
Canova, Fabio, and Luca Sala. 2009. “Back to square one: Identification issues in DSGE models.” Journal of Monetary Economics 56 (4): 431–49. https://doi.org/10.1016/j.jmoneco.2009.03.014.
Carrella, Ernesto, Richard M. Bailey, and Jens Koed Madsen. 2019. “Repeated discrete choices in geographical agent based models with an application to fisheries.” Environmental Modelling and Software 111 (January): 204–30. https://doi.org/10.1016/j.envsoft.2018.08.023.
Carruthers, Thomas R., André E. Punt, Carl J. Walters, Alec MacCall, Murdoch K. McAllister, Edward J. Dick, and Jason Cope. 2014. “Evaluating methods for setting catch limits in data-limited fisheries.” Fisheries Research 153 (May): 48–68. https://doi.org/10.1016/j.fishres.2013.12.014.
Cramer, Simon, and Torsten Trimborn. 2019. “Stylized Facts and Agent-Based Modeling.” https://stooq.com.
Franke, Reiner, and Frank Westerhoff. 2016. “Why a simple herding model may generate the stylized facts of daily returns: explanation and estimation.” Journal of Economic Interaction and Coordination 11 (1): 1–34. https://doi.org/10.1007/s11403-014-0140-6.
Froese, Rainer, Henning Winker, Gianpaolo Coro, Nazli Demirel, Athanassios C Tsikliras, Donna Dimarchopoulou, Giuseppe Scarcella, Wolfgang Nikolaus Probst, Manuel Dureuil, and Daniel Pauly. 2019a. “On the pile-up effect and priors for Linf and M/K: response to a comment by Hordyk et al. on ‘A new approach for estimating stock status from length frequency data’.” ICES Journal of Marine Science 76 (2): 461–65. https://doi.org/10.1093/icesjms/fsy199.
———. 2019b. “A new approach for estimating stock status from length frequency data.” ICES Journal of Marine Science 76 (1): 350–51. https://doi.org/10.1093/icesjms/fsy139.
Goudriaan, J. 1986. “Boxcartrain Methods for Modelling of Ageing, Development, Delays and Dispersion.” In The Dynamics of Physiologically Structured Populations, edited by J. A. J. Metz and O. Dierkmann, 1st ed., 453–73. Berlin: Springer.
Grazzini, Jakob, Matteo G. Richiardi, and Mike Tsionas. 2017. “Bayesian estimation of agent-based models.” Journal of Economic Dynamics and Control 77 (April): 26–47. https://doi.org/10.1016/j.jedc.2017.01.014.
Grimm, Volker, Eloy Revilla, Uta Berger, Florian Jeltsch, Wolf M. Mooij, Steven F. Railsback, Hans Hermann Thulke, Jacob Weiner, Thorsten Wiegand, and Donald L. DeAngelis. 2005. “Pattern-oriented modeling of agent-based complex systems: Lessons from ecology.” https://doi.org/10.1126/science.1116681.
Hartig, F., J. M. Calabrese, B. Reineking, T. Wiegand, and A. Huth. 2011. “Statistical inference for stochastic simulation models - theory and application.” Ecology Letters 14 (8): 816–27. https://doi.org/10.1111/j.1461-0248.2011.01640.x.
Hordyk, Adrian, Kotaro Ono, Keith Sainsbury, Neil Loneragan, and Jeremy Prince. 2014. “Some explorations of the life history ratios to describe length composition, spawning-per-recruit, and the spawning potential ratio.” In ICES Journal of Marine Science, 72:204–16. 1. https://doi.org/10.1093/icesjms/fst235.
Hordyk, Adrian R., and Thomas R. Carruthers. 2018. “A quantitative evaluation of a qualitative risk assessment framework: Examining the assumptions and predictions of the Productivity Susceptibility Analysis (PSA).” PLoS ONE 13 (6). https://doi.org/10.1371/journal.pone.0198298.
Hordyk, Adrian R, Jeremy D Prince, Thomas R Carruthers, and Carl J Walters. 2019. “Comment on "a new approach for estimating stock status from length frequency data" by Froese et al. (2018).” ICES Journal of Marine Science 76 (2): 457–60. https://doi.org/10.1093/icesjms/fsy168.
Makino, Mitsutaku, Shingo Watari, Taro Hirose, Kentaro Oda, Masahito Hirota, Atsushi Takei, Michio Ogawa, and Hiroshi Horikawa. 2017. “A transdisciplinary research of coastal fisheries co-management: the case of the hairtail Trichiurus japonicus trolling line fishery around the Bungo Channel, Japan.” Fisheries Science 83 (6): 853–64. https://doi.org/10.1007/s12562-017-1141-x.
Mangel, Marc, Alec D. MacCall, Jon Brodziak, E. J. Dick, Robyn E. Forrest, Roxanna Pourzand, and Stephen Ralston. 2013. “A perspective on steepness, reference points, and stock assessment.” Edited by Kenneth Rose. Canadian Journal of Fisheries and Aquatic Sciences 70 (6): 930–40. https://doi.org/10.1139/cjfas-2012-0372.
Manski, Charles F. 2003. Partial identification of probability distributions. New York, New York, USA: Springer.
———. 2007. Identification for Prediction and Decision. https://books.google.com.tw/books/about/Identification{\_}for{\_}Prediction{\_}and{\_}Decisi.html?id=KcW-820XgeQC{\&}printsec=frontcover{\&}source=kp{\_}read{\_}button{\&}redir{\_}esc=y{\#}v=onepage{\&}q{\&}f=false.
Martell, Steven, and Rainer Froese. 2013. “A simple method for estimating MSY from catch and resilience.” Fish and Fisheries 14 (4): 504–14. https://doi.org/10.1111/j.1467-2979.2012.00485.x.
Mildenberger, Tobias Karl, Marc Hollis Taylor, and Matthias Wolff. 2017. “<b>TropFishR</b> : an R package for fisheries analysis with length-frequency data.” Edited by Samantha Price. Methods in Ecology and Evolution 8 (11): 1520–7. https://doi.org/10.1111/2041-210X.12791.
Pearl, Judea. 2014. “The Causal Foundations of Structural Equation Modeling.” Handbook of Structural Equation Modeling, no. June: 68–91. http://ftp.cs.ucla.edu/pub/stat{\_}ser/r370.pdf.
———. 2017. “A Linear ‘Microscope’ for Interventions and Counterfactuals DOI.” J. Causal Infer. https://doi.org/10.1515/jci-2017-0003.
Quinn, Terrance J., and Richard B. Deriso. 1999. Quantitative Fish Dynamics. Oxford University Press.
Russell, Stuart, and Peter Norvig. 2009. Artificial Intelligence: a modern approach. 3rd ed. Upper Saddle River, NJ, USA: Prentice Hall Press.
Seijo, Juan Carlos, Omar Defeo, and Silvia Salas. 1998. Fisheries bioeconomics Theory, modelling and management. FAO. http://www.fao.org/3/W6914E/W6914E02.htm.
Sisson, Scott A. 2018. Handbook of Approximate Bayesian Computation. https://doi.org/10.1201/9781315117195.
Smith, Vernon L. 1969. “On Models of Commercial Fishing.” Journal of Political Economy 77 (2): 181–98. https://doi.org/10.1086/259507.
Sparre, Per, and Siebren C. Venema. 1998. Introduction to tropical fish stock assessment Introduction to tropical fish stock assessment Part I : Manual. Vol. 306. https://ci.nii.ac.jp/naid/10030499433/.
Thorson, James T., Stephan B. Munch, Jason M. Cope, and Jin Gao. 2017. “Predicting life history parameters for all fishes worldwide.” 8. Vol. 27. https://doi.org/10.1002/eap.1606.
Varian, Hal. 1989. “What Use is Economic Theory?”
Watari, Shingo, Syunji Tokumitsu, Taro Hirose, Michio Ogawa, and Mitsutaku Makino. 2017. “Stock structure and resource management of hairtail Trichiurus japonicus based on seasonal broods around the Bungo Channel, Japan.” Fisheries Science 83 (6): 865–78. https://doi.org/10.1007/s12562-017-1140-y.
Zhang, Jingjing, Todd E. Dennis, Todd J. Landers, Elizabeth Bell, and George L. W. Perry. 2017. “Linking individual-based and statistical inferential models in movement ecology: A case study with black petrels (Procellaria parkinsoni).” Ecological Modelling 360: 425–36. https://doi.org/10.1016/j.ecolmodel.2017.07.017.
This takes approximately 30 minutes on an average 2018 laptop↩︎